/**************************************************************************
* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
*                                                                        *
* Author: The ALICE Off-line Project.                                    *
* Contributors are mentioned in the code where appropriate.              *
*                                                                        *
* Permission to use, copy, modify and distribute this software and its   *
* documentation strictly for non-commercial purposes is hereby granted   *
* without fee, provided that the above copyright notice appears in all   *
* copies and that both the copyright notice and this permission notice   *
* appear in the supporting documentation. The authors make no claims     *
* about the suitability of this software for any purpose. It is          *
* provided "as is" without express or implied warranty.                  *
**************************************************************************/

// $Id: AliAODMuonReplicator0.cxx 56492 2012-05-15 18:42:47Z pcrochet $

//
// Implementation of a branch replicator 
// to produce slim muon and dimuon aods.
//
// This replicator is in charge of replicating the tracks,vertices,dimuons
// (and vzero and tzero) branches of the standard AOD into muon AODs 
// (AliAOD.Muons.root and AliAOD.Dimuons.root)
// 
// The tracks are filtered so that only muon tracks (and only muon tracks
// that pass the trackCut if present) make it to the output aods
//
// The vertices are filtered so that only the primary vertices make it
// to the output aods.
//
// The dimuons are recreated here, according to the set of tracks
// that pass the trackCut (that set may be the same as the input set,
// but to be 100% safe, we always recreate the dimuons).
// 
// Author: L. Aphecetche (Subatech)


class AliESDv0;
class AliESDVertex;
class AliAODVertex;
//class AliAOD3LH;

#include "AliAODDimuon.h"
#include "AliAODEvent.h"
#include "AliAODMCHeader.h"
#include "AliAODMCParticle.h"
#include "AliAODTZERO.h"
#include "AliAODTrack.h"
#include "AliAODVZERO.h"
#include "AliAnalysisCuts.h"
#include "TF1.h"
#include "AliExternalTrackParam.h"
#include "AliESDv0.h"
#include "AliAODv0.h"
#include "AliPIDResponse.h"
//#include "AliAOD3LH.h"
#include <iostream>
#include <cassert>
#include "AliESDtrack.h"
#include "AliAODMuonReplicator0.h"

using std::cout;
using std::endl;

ClassImp(AliAODMuonReplicator0)

//_____________________________________________________________________________
AliAODMuonReplicator0::AliAODMuonReplicator0(const char* name, const char* title,
					     AliAnalysisCuts* trackCut,
					     AliAnalysisCuts* vertexCut,
					     Int_t mcMode, 
					     Int_t nsigmaTrk1, Int_t partType1,
					     Int_t nsigmaTrk2, Int_t partType2)
:AliAODBranchReplicator(name,title), 
  fTrackCut(trackCut), fTracks(0x0), 
  fVertexCut(vertexCut), fVertices(0x0), 
  fDimuons(0x0),
  fVZERO(0x0),
  fTZERO(0x0),
  fList(0x0),
  fMCParticles(0x0),
  fMCHeader(0x0),
  fnSigmaTrk1(nsigmaTrk1),
  fnSigmaTrk2(nsigmaTrk2),
  fpartType1(partType1),
  fpartType2(partType2),
  fMCMode(mcMode),
  
  fLabelMap(),
  fParticleSelected(),
  fAODTrackRefs(0x0)
{
  // default ctor
}

//_____________________________________________________________________________
AliAODMuonReplicator0::~AliAODMuonReplicator0()
{
  // dtor
  delete fTrackCut;
  delete fVertexCut;
  delete fList;
}

//_____________________________________________________________________________
void AliAODMuonReplicator0::SelectParticle(Int_t i)
{
  // taking the absolute values here, need to take care 
  // of negative daughter and mother
  // IDs when setting!
  
  if (!IsParticleSelected(TMath::Abs(i)))
  {
    fParticleSelected.Add(TMath::Abs(i),1);    
  }
}

//_____________________________________________________________________________
Bool_t AliAODMuonReplicator0::IsParticleSelected(Int_t i)  
{
  // taking the absolute values here, need to take 
  // care with negative daughter and mother
  // IDs when setting!
  return (fParticleSelected.GetValue(TMath::Abs(i))==1);
}


//_____________________________________________________________________________
void AliAODMuonReplicator0::CreateLabelMap(const AliAODEvent& source)
{  
  //
  // this should be called once all selections are done 
  //
  
  fLabelMap.Delete();
  
  TClonesArray* mcParticles = static_cast<TClonesArray*>(source.FindListObject(AliAODMCParticle::StdBranchName()));
  
  Int_t i(0);
  Int_t j(0);
  
  TIter next(mcParticles);
  
  while ( next() ) 
  {
    if (IsParticleSelected(i))
    {
      fLabelMap.Add(i,j++);
    }
    ++i;
  }
}

//_____________________________________________________________________________
Int_t AliAODMuonReplicator0::GetNewLabel(Int_t i) 
{
  // Gets the label from the new created Map
  // Call CreatLabelMap before
  // otherwise only 0 returned
  return fLabelMap.GetValue(TMath::Abs(i));
}

//_____________________________________________________________________________
void AliAODMuonReplicator0::FilterMC(const AliAODEvent& source)
{
  // Filter MC information

  fMCHeader->Reset();
  fMCParticles->Clear("C");

  AliAODMCHeader* mcHeader(0x0);
  TClonesArray* mcParticles(0x0);
  
  fParticleSelected.Delete();
  
  if ( fMCMode>=2 && !fTracks->GetEntries() ) return;
  // for fMCMode==1 we only copy MC information for events where there's at least one muon track
    
  mcHeader = static_cast<AliAODMCHeader*>(source.FindListObject(AliAODMCHeader::StdBranchName()));
  
  if ( mcHeader ) 
    {
      *fMCHeader = *mcHeader;
    }
  
  mcParticles = static_cast<TClonesArray*>(source.FindListObject(AliAODMCParticle::StdBranchName()));
  
  if ( mcParticles && fMCMode>=2 )
    {
      // loop on (kept) muon tracks to find their ancestors
      TIter nextMT(fTracks);
      AliAODTrack* mt;
    
      while ( ( mt = static_cast<AliAODTrack*>(nextMT()) ) )
	{
	  Int_t label = mt->GetLabel();
      
	  while ( label >= 0 ) 
	    {
	      SelectParticle(label);
	      AliAODMCParticle* mother = static_cast<AliAODMCParticle*>(mcParticles->UncheckedAt(label));
	      if (!mother)
		{
		  AliError("Got a null mother ! Check that !");
		  label = -1;
		}
	      else
		{
		  label = mother->GetMother();
		}
	    }
	}
    
      CreateLabelMap(source);
    
      // Actual filtering and label remapping (shamelessly taken for the implementation of AliAODHandler::StoreMCParticles)
      TIter nextMC(mcParticles);
      AliAODMCParticle* p;
      Int_t nmc(0);
      Int_t nmcout(0);
    
      while ( ( p = static_cast<AliAODMCParticle*>(nextMC()) ) )
	{
	  AliAODMCParticle c(*p);
      
	  if ( IsParticleSelected(nmc) )
	    {
	      // 
	      Int_t d0 =  p->GetDaughter(0);
	      Int_t d1 =  p->GetDaughter(1);
	      Int_t m =   p->GetMother();
        
	      // other than for the track labels, negative values mean
	      // no daughter/mother so preserve it
        
	      if(d0<0 && d1<0)
		{
		  // no first daughter -> no second daughter
		  // nothing to be done
		  // second condition not needed just for sanity check at the end
		  c.SetDaughter(0,d0);
		  c.SetDaughter(1,d1);
		} else if(d1 < 0 && d0 >= 0) 
		{
		  // Only one daughter
		  // second condition not needed just for sanity check at the end
		  if(IsParticleSelected(d0))
		    {
		      c.SetDaughter(0,GetNewLabel(d0));
		    } else 
		    {
		      c.SetDaughter(0,-1);
		    }
		  c.SetDaughter(1,d1);
		}
	      else if (d0 > 0 && d1 > 0 )
		{
		  // we have two or more daughters loop on the stack to see if they are
		  // selected
		  Int_t d0tmp = -1;
		  Int_t d1tmp = -1;
		  for (int id = d0; id<=d1;++id)
		    {
		      if (IsParticleSelected(id))
			{
			  if(d0tmp==-1)
			    {
			      // first time
			      d0tmp = GetNewLabel(id);
			      d1tmp = d0tmp; // this is to have the same schema as on the stack i.e. with one daugther d0 and d1 are the same 
			    }
			  else d1tmp = GetNewLabel(id);
			}
		    }
		  c.SetDaughter(0,d0tmp);
		  c.SetDaughter(1,d1tmp);
		} else 
		{
		  AliError(Form("Unxpected indices %d %d",d0,d1));
		}
        
	      if ( m < 0 )
		{
		  c.SetMother(m);
		} else 
		{
		  if (IsParticleSelected(m)) 
		    {
		      c.SetMother(GetNewLabel(m));              
		    }
		  else 
		    {
		      AliError(Form("PROBLEM Mother not selected %d", m));              
		    }
		}
        
	      new ((*fMCParticles)[nmcout++]) AliAODMCParticle(c);
	    }
      
	  ++nmc;        
	}      
    
      // now remap the tracks...
    
      TIter nextTrack(fTracks);
      AliAODTrack* t;
    
      while ( ( t = static_cast<AliAODTrack*>(nextTrack()) ) )
	{
	  t->SetLabel(GetNewLabel(t->GetLabel()));
	}
    
    }
  else if ( mcParticles ) 
    {
      // simple copy of input MC particles to ouput MC particles
      TIter nextMC(mcParticles);
      AliAODMCParticle* p;
      Int_t nmcout(0);
    
      while ( ( p = static_cast<AliAODMCParticle*>(nextMC()) ) )
	{
	  new ((*fMCParticles)[nmcout++]) AliAODMCParticle(*p);
	}
    }
  
  AliDebug(1,Form("input mc %d output mc %d",
                  mcParticles ? mcParticles->GetEntries() : 0,
                  fMCParticles ? fMCParticles->GetEntries() : 0));
  
}

//_____________________________________________________________________________
TList* AliAODMuonReplicator0::GetList() const
{
  // return (and build if not already done) our internal list of managed objects
  if (!fList)
    {
      // fTracks = new TClonesArray("AliAODTrack",30);
      // fTracks->SetName("tracks");    

      fTracks = new TClonesArray("AliAODv0",30);
      // fTracks->SetName("SecondaryVertices");    
      
      //  fTracks = new TClonesArray("AliESDv0",30);
      fTracks->SetName("SecondaryVertices");    

    
      fVertices = new TClonesArray("AliAODVertex",2);
      fVertices->SetName("vertices");    
    
      //    fDimuons = new TClonesArray("AliAODDimuon",30);
      fDimuons = new TClonesArray("AliAODTrack",30);
      fDimuons->SetName("Nuclei");
    
      fVZERO = new TClonesArray("AliAODTrack",30);
      fVZERO->SetName("DaughterTracks");

      fTZERO = new AliAODTZERO;
    
      fList = new TList;
      fList->SetOwner(kTRUE);
    
      fList->Add(fTracks);
      fList->Add(fVertices);
      fList->Add(fDimuons);
      fList->Add(fVZERO);
      
      //fList->Add(fTZERO);
    
      if ( fMCMode > 0 )
	{
	  fMCHeader = new AliAODMCHeader;    
	  fMCParticles = new TClonesArray("AliAODMCParticle",1000);
	  fMCParticles->SetName(AliAODMCParticle::StdBranchName());
	  fList->Add(fMCHeader);
	  fList->Add(fMCParticles);
	}
    }
  return fList;
}

//_____________________________________________________________________________
void AliAODMuonReplicator0::ReplicateAndFilter(const AliAODEvent& source)
{
  // Replicate (and filter if filters are there) the relevant parts we're interested in AODEvent
  
  assert(fTracks!=0x0);
  
  *fTZERO = *(source.GetTZEROData());

  fTracks->Clear("C");			
    
  fDimuons->Clear("C");

  fVZERO->Clear("C");
  
  TIter next(source.GetTracks());
  AliAODTrack* t;
  Int_t ntracks(0);
  Int_t ndimuons(0);

  Int_t input(0);

 
  AliAODVertex *vtx = source.GetPrimaryVertex();
  
  if(vtx->GetNContributors()<1) {
    
    // SPD vertex cut
    vtx =source.GetPrimaryVertexSPD();
    
    if(vtx->GetNContributors()<1) {
      Info("AliAnalysisTaskHelium3Pi","No good vertex, skip event");
      return; // NO GOOD VERTEX, SKIP EVENT 
    }
  }

  // Double_t xPrimaryVertex=0.,yPrimaryVertex=0.,zPrimaryVertex=0.;
  // xPrimaryVertex=vtx->GetX();
  // yPrimaryVertex=vtx->GetY();
  // zPrimaryVertex=vtx->GetZ();  

  Double_t lMagneticField=source.GetMagneticField();
  
  Double_t position[3];
  for(Int_t i=0;i<3;i++)position[i]=0;
  
  Float_t impactXY=-999, impactZ=-999;

  printf("Sto eseguendo il muon Replicator...\n");

  TF1 *foPion=new TF1("foPion", "[0]*([1]*TMath::Power(TMath::Sqrt(1 + (x/0.13957)*(x/0.13957))/(x/0.13957) , [3]) - 1 - TMath::Power(TMath::Sqrt(1 + (x/0.13957)*(x/0.13957))/(x/0.13957) , [3])*TMath::Log([2] + 1/TMath::Power((x/0.13957), [4])))",0.01,20);
  foPion->SetParameters(4.1,8.98482806165147636e+00,1.54000000000000005e-05,2.30445734159456084e+00,2.25624744086878559e+00);
  
  Double_t TrackNumber = source.GetNumberOfTracks();
  Int_t    countPionBuoni=0;
  
  //Tracks arrays

  TArrayI Pions(TrackNumber);     //Pions                                                                          

  Int_t nPions=0;

  TArrayI He(TrackNumber);        //Helium3
  Int_t nHe=0;
  
  fAODTrackRefs = new TRefArray(TrackNumber);
 
  ULong_t  status;
  Bool_t   isTPC, isTOF;

  for(Int_t j=0; j<TrackNumber; j++){

    AliAODTrack *aodtrack=source.GetTrack(j);
    
    //----------------------------------------
    // To be implemented into a method IsTrackAccepted 
    
    status  = (ULong_t)aodtrack->GetStatus();
    // isTPC   = ((status & AliAODTrack::kTPCin)  != 0);
    // isTOF   = (((status & AliAODTrack::kTOFout) != 0) && ((status & AliAODTrack::kTIME) != 0));
    
    // //     if((status&AliESDtrack::kITSrefit)==0) continue;// Require itsrefit
    if((status&AliESDtrack::kTPCrefit)==0) continue;// Require tpcrefit
    
    //    if(aodtrack->GetTPCNcls() < 60) continue; // TPC Ncluster cut to be fixed
    if(aodtrack->GetTPCNcls() < 80) continue; // TPC Ncluster cut to be fixed
    if(aodtrack->Chi2perNDF() > 5) continue;

    //-----------------------------------------

    Double_t mom = aodtrack->GetDetPid()->GetTPCmomentum();
    // Double_t mom =  aodtrack->P();
    
    //if(mom<0.1)continue;
    if(mom<0.150)continue;
    
    // Double_t nSigmaNegPion = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(aodtrack,(AliPID::EParticleType)fpartType1)); //to be changed

    Double_t nSigmaNegPion = TMath::Abs((aodtrack->GetTPCsignal() - foPion->Eval(mom))/foPion->Eval(mom))/0.07;
    
    if ( (nSigmaNegPion <= fnSigmaTrk1)){
      Pions[nPions++]=j;
    }
    
    // Double_t nSigmaNuclei = fPIDResponse->NumberOfSigmasTPC(aodtrack,(AliPID::EParticleType)fpartType2);  //To be changed
    // if(nSigmaNuclei>-fnSigmaTrk2 && aodtrack->GetTPCsignal()<1000 && mom>0.2){ 

    Double_t triggerDeDx = 4*AliExternalTrackParam::BetheBlochAleph((mom*2)/(0.938*3),1.0288,31.9806,5.04114e-11,2.13096,2.38541);
    
    if(aodtrack->GetTPCsignal() > triggerDeDx && aodtrack->GetTPCsignal()<5000 && mom>0.2){
      He[nHe++]=j;
      
      new((*fDimuons)[ndimuons++]) AliAODTrack(*aodtrack);
      //cout<<"He tracks : "<<aodtrack->GetTPCNcls()<<endl;
    }
  } 

  Pions.Set(nPions);
  He.Set(nHe);

  // cout<<"Numero di pi: "<<nPions<<"; N di 3He: "<<nHe<<" Dimensione pi vector"<<Pions.GetSize()<<" and 3He: "<<He.GetSize()<<endl;  
  
  //------------------------------
  
  ReconstructAndSaveSecondaryVertex(source,Pions,He);
  cout<<"  fAODTrackRefs " <<fAODTrackRefs<<endl;

  //------------------------------

  // to be understand : do we need the information of the daughters? if yes we sould do someting like
  // $ALICE_ROOT/ANALYSIS/AliAnalysisTaskESDfilter.cxx :: Convert v0(AliESDevent) or not??
  // forse si se vogliamo poi fare dei controlli aggiuntivi sulle tracce
  // NB pero le abbiamo gia convertite, quindi sono solo piu da aggiungere in qualche modo
    
  //---------------------------------------------------------------------------------
  
  assert(fVertices!=0x0);
  fVertices->Clear("C");
  TIter nextV(source.GetVertices());
  AliAODVertex* v;
  Int_t nvertices(0);
  
  while ( ( v = static_cast<AliAODVertex*>(nextV()) ) )
    {
      if ( !fVertexCut || fVertexCut->IsSelected(v) ) 
	{
	  AliAODVertex* tmp = v->CloneWithoutRefs();
	  AliAODVertex* copiedVertex = new((*fVertices)[nvertices++]) AliAODVertex(*tmp);
	  
	  // to insure the main vertex retains the ncontributors information
	  // (which is otherwise computed dynamically from
	  // references to tracks, which we do not keep in muon aods...)
	  // we set it here
	  
	  copiedVertex->SetNContributors(v->GetNContributors()); 
	  delete tmp;
	}
    }
  
  printf("....Ho eseguito il muon Replicator...\n");
  
  AliDebug(1,Form("input mu tracks=%d tracks=%d vertices=%d ndimuons=%d",
                  input,fTracks->GetEntries(),fVertices->GetEntries(),fDimuons->GetEntries())); 
  
  cout<<"....Sotto debug..."<<endl;
  
  // Finally, deal with MC information, if needed
  
  if ( fMCMode > 0 )
    {
      FilterMC(source);      
    }

  delete fAODTrackRefs; fAODTrackRefs=0x0;

  cout<<"Fine Replicator"<<endl;
}

//___________________________________________________________________________________________

//void AliAODMuonReplicator0::ReconstructAndSaveSecondaryVertex(const AliAODEvent& source, TArrayI &Pions, TArrayI &He)

// {

//   Int_t nV0(0);
//   Int_t ntracksD(0);

//   AliAODVertex *vtx = source.GetPrimaryVertex();
//   AliESDv0 *vertex ;

//   //! A set of very loose parameters for cuts 

//   Double_t fgChi2max=33.;     //! max chi2
//   Double_t fgDNmin=0.05;      //! min imp parameter for the 1st daughter = 500um
//   Double_t fgDPmin=0.05;      //! min imp parameter for the 2nd daughter = 500um
//   Double_t fgDCAmax=1.;       //! max DCA between the daughter tracks in cm
//   Double_t fgCPAmin=0.9;      //! min cosine of V0's pointing angle
//   Double_t fgRmin=0.2;        //! min radius of the fiducial volume //original
//   //  Double_t fgRmin=0.1;        //! min radius of the fiducial volume = 1 mm 
//   Double_t fgRmax=200.;       //! max radius of the fiducial volume = 2 m

//   //------------------------------------------
  
//   Double_t dcaPionToPV[2] = {1e6,1e6};
//   Double_t covarPionToPV[3] = {1e6,1e6,1e6};
    
//   Double_t dcaHeToPV[2] = {1e6,1e6};
//   Double_t covarHeToPV[3] = {1e6,1e6,1e6};

//   Double_t DcaPionToPrimVertex=0;
//   Double_t DcaHeToPrimVertex=0;

//   Double_t pos[3] = { 0. };      
//   Double_t chi2(0.0);
//   Double_t covVtx[6] = { 0. };

//   //---------------------------------------------------------------------------------------------

//   AliAODTrack  *PionTrk = 0x0;
//   AliAODTrack  *HeTrack = 0x0;
  
//   AliAODv0 *myvertex;

//   Int_t nHe    = He.GetSize();
//   Int_t nPions = Pions.GetSize();

//   Int_t posFromV0 = 0;
//   Int_t negFromV0 = 0;

//   cout<<"nHe: "<<nHe<<" npions: "<<nPions<<endl;

//   for (Int_t i=0; i<nHe; i++){                            //! He Tracks Loop

//     Int_t He3idx=He[i];
    
//     HeTrack=source.GetTrack(He3idx);
     
//     if (HeTrack){
//       HeTrack->PropagateToDCA(source.GetPrimaryVertex(),source.GetMagneticField(),100.,dcaHeToPV,covarHeToPV);
//       DcaHeToPrimVertex = TMath::Sqrt(dcaHeToPV[0]*dcaHeToPV[0]+dcaHeToPV[1]*dcaHeToPV[1]);
//     }
    
//     //-----

//     for (Int_t k=0; k < nPions; k++) {                           //! Pion Tracks Loop
      
//       Int_t PionIdx=Pions[k];
      
//       PionTrk=source.GetTrack(PionIdx);
      
//       if (PionTrk) {
// 	PionTrk->PropagateToDCA(source.GetPrimaryVertex(),source.GetMagneticField(),100.,dcaPionToPV,covarPionToPV);
// 	DcaPionToPrimVertex = TMath::Sqrt(dcaPionToPV[0]*dcaPionToPV[0]+dcaPionToPV[1]*dcaPionToPV[1]);	
//       }
      
//       //   from   $ALICE_ROOT/STERR/ESD/AliV0vertexer

//       if ( DcaPionToPrimVertex < fgDNmin)                
// 	if ( DcaHeToPrimVertex < fgDNmin) continue;      
      
//       Double_t xn, xp;
//       Double_t dca=0.;
      
//       AliExternalTrackParam trackInHe3TPC,trackInNegPion;

//       trackInHe3TPC.CopyFromVTrack(HeTrack);  //copy constructor for AliExternalTrackParam object
//       trackInNegPion.CopyFromVTrack(PionTrk);

//       if(HeTrack==PionTrk)continue;
//       if(PionIdx==He3idx)continue;

//       // AliESDtrack *esdtrack1 = new AliESDtrack((AliVTrack*)HeTrack);
//       // AliESDtrack *esdtrack2 = new AliESDtrack((AliVTrack*)PionTrk);

//       // if(esdtrack1==esdtrack2)continue;
//       // if(!esdtrack1->GetCovariance())continue; 
//       // if(!esdtrack2->GetCovariance())continue; 
      
//       // if(!esdtrack1->GetSigmaY2())continue;
//       // if(!esdtrack1->GetSigmaZ2())continue;
//       // if(!esdtrack2->GetSigmaY2())continue;
//       // if(!esdtrack2->GetSigmaZ2())continue;

//       //cout<<trackInHe3TPC.CopyFromVTrack(static_cast<AliAODTrack*>(HeTrack))<<endl;
      
//       dca = trackInHe3TPC.GetDCA(&trackInNegPion,source.GetMagneticField(),xn,xp);  //!dca tracks
//       cout<<"DCA: "<<dca<<endl;
	 
//       //      dca = 0.5;
      
    
//       //    if(! trackInHe3TPC.CopyFromVTrack(static_cast<AliAODTrack*>(HeTrack)))continue;
//       // cout<<trackInHe3TPC<<endl;

//       //      if(trackInHe3TPC.GetDCA(&trackInNegPion,source.GetMagneticField(),xn,xp)<=0)continue;
      
//       //      Double_t dca = trackInHe3TPC.GetDCA(&trackInNegPion,source.GetMagneticField(),xn,xp);  //!dca tracks

//       //      if(esdtrack1->GetDCA(esdtrack2,source.GetMagneticField(),xn,xp)<=0)continue;
  
//       //      Double_t dca = esdtrack1->GetDCA(esdtrack2,source.GetMagneticField(),xn,xp);
//       cout<<dca<<endl;
//       cout<<DcaPionToPrimVertex<<endl;
//       cout<<DcaHeToPrimVertex<<endl;
//       // cout<<esdtrack1->GetD(vtx->GetX(),vtx->GetY(),source.GetMagneticField())<<endl;
//       // cout<<esdtrack2->GetD(vtx->GetX(),vtx->GetY(),source.GetMagneticField())<<endl;
      
//       //      if(!dca)continue;
  
//       if ((xn+xp) > 2*fgRmax) continue;
//       if ((xn+xp) < 2*fgRmin) continue;
      
//       if (dca > fgDCAmax) continue;
//       // if (dca <= 0) continue;
      
       
//       //From AliV0Vertex
      
//       Bool_t corrected=kFALSE;
//       if ((trackInNegPion.GetX() > 3.) && (xn < 3.)) {
// 	//correct for the beam pipe material
// 	corrected=kTRUE;
//       }
//       if ((trackInHe3TPC.GetX() > 3.) && (xp < 3.)) {
// 	//correct for the beam pipe material
// 	corrected=kTRUE;
//       }
      
//       // if (corrected) {
//       // 	dca=esdtrack1->GetDCA(esdtrack2,source.GetMagneticField(),xn,xp);//trackInNegPion.GetDCA(&trackInHe3TPC,source.GetMagneticField(),xn,xp);
//       // 	if (dca > fgDCAmax) continue;
//       // 	if ((xn+xp) > 2*fgRmax) continue;
//       // 	if ((xn+xp) < 2*fgRmin) continue;
//       // }
      

//       //===================//
//       // "V0" with Tracks  //
//       //===================//

//       trackInNegPion.PropagateTo(xn,source.GetMagneticField()); 
//       trackInHe3TPC.PropagateTo(xp,source.GetMagneticField());
      
//       vertex = new AliESDv0(trackInNegPion,PionIdx,trackInHe3TPC,He3idx);
//       vertex->SetDcaV0Daughters(dca);
//       vertex->SetV0CosineOfPointingAngle(vertex->GetV0CosineOfPointingAngle(vtx->GetX(),vtx->GetY(),vtx->GetZ()));
//       vertex->GetXYZ(pos[0],pos[1],pos[2]);
   
//       chi2 = vertex->GetChi2V0(); // = chi2/NDF since NDF = 2*2-3
//       vertex->GetPosCov(covVtx);

//       if(vertex->GetV0CosineOfPointingAngle(vtx->GetX(),vtx->GetY(),vtx->GetZ())<fgCPAmin)continue;

//       Double_t rDca3LHDaughters = dca;
//       Double_t rDca3LHToPrimVertex = vertex->GetD(vtx->GetX(),vtx->GetY(),vtx->GetZ());
  
//       Double_t rMomPos[3]={0.,0.,0.};
//       Double_t rMomNeg[3]={0.,0.,0.};
      
//       Double_t rDCA[2]={0.,0.}; //to be fixed!
  
//       //ramona 
      
//       rDCA[0] =  DcaHeToPrimVertex;
//       rDCA[1] =  DcaPionToPrimVertex;
      
//       //to be understand if pos is 3He and neg is pi in each case
  
//       vertex->GetPPxPyPz(rMomPos[0], rMomPos[1], rMomPos[2]); //is momentum at SV!!
//       vertex->GetNPxPyPz(rMomNeg[0], rMomNeg[1], rMomNeg[2]);
  
//       const Double_t rMomPos1[3]={2*rMomPos[0], 2*rMomPos[1], 2*rMomPos[2]};
//       const Double_t rMomNeg1[3]={rMomNeg[0], rMomNeg[1], rMomNeg[2]};

//       // myvertex=new AliAODv0(vtx,rDca3LHDaughters,rDca3LHToPrimVertex,rMomPos1,rMomNeg1,rDCA);
//       // new((*fTracks)[nV0++]) AliAODv0(*myvertex);

//       // cout<<"cpa: "<<vertex->GetV0CosineOfPointingAngle(vtx->GetX(),vtx->GetY(),vtx->GetZ())<<endl;
//       // cout<<"momPospx: "<<rMomPos[0]<<" 1 "<<rMomPos1[0]<<endl;
//       // cout<<"momNegpx: "<<rMomNeg[0]<<" 1 "<<rMomNeg1[0]<<endl;

//       // //-----

//       //Conversion from v0 Coneversion
      
//       posFromV0 = vertex->GetPindex();
//       negFromV0 = vertex->GetNindex();

     
//       //Additional infos
      
//       AliAODVertex * vV0 =  new AliAODVertex(pos,
//       					     covVtx,
//       					     chi2,
//       					     vtx,
//       					     nV0,
//       					     AliAODVertex::kV0);
    
//       // esdV0Pos->GetPxPyPz(momPos);
//       // esdV0Pos->GetXYZ(pos);
//       // esdV0Pos->GetCovarianceXYZPxPyPz(covTr);
//       // esdV0Pos->GetESDpid(pid);
//       // esdV0Pos->GetTOFLabel(tofLabel);
      
      
//       //      HeTrack = static_cast<AliAODTrack*>(fAODTrackRefs->At(posFromV0));
      
      
//       // cout<<"\n\n---------> Sec vtx: "<<vertex<<endl;
//       // cout<<"           Pos idx from ESD v0: "<<posFromV0<<" He Tracks: "<<HeTrack<<endl;
//       // cout<<"           Neg idx from ESD v0: "<<negFromV0<<" Pi Tracks: "<<PionTrk<<endl<<endl;
     
//       fAODTrackRefs->AddAt(HeTrack,posFromV0);
//       fAODTrackRefs->AddAt(PionTrk,negFromV0);
     
  
//       // vertex->SetIndex(0,posFromV0);
//       // vertex->SetIndex(1,negFromV0);


//       vV0->AddDaughter(HeTrack);
//       vV0->AddDaughter(PionTrk);
//       vV0->SetPosition(pos[0],pos[1],pos[2]);
      
//       // cout<<"           Figlie from vV0: "<<vV0->GetNDaughters()<<endl;
//       // cout<<"           Pos idx from vV0: "<<vV0->GetDaughter(0)<<endl;
//       // cout<<"           Neg idx from vV0: "<<vV0->GetDaughter(1)<<endl<<endl;

      
//       //era questo
//       myvertex=new AliAODv0(vV0,rDca3LHDaughters,rDca3LHToPrimVertex,rMomPos1,rMomNeg1,rDCA);
//       new((*fTracks)[nV0++]) AliAODv0(*myvertex);

//       //new((*fTracks)[nV0++]) AliESDv0(*vertex);

      
//       AliAODv0* myvertexv0 =  new AliAODv0(vV0, 
// 					   rDca3LHDaughters,
// 					   rDca3LHToPrimVertex,
// 					   rMomPos1,
// 					   rMomNeg1,
// 					   rDCA);

//       // cout<<"           Pos idx from AODV0 v0:  "<<myvertexv0->GetPosID()<<endl;
//       // cout<<"           Neg idx from AODV0 v0:  "<<myvertexv0->GetNegID()<<endl;
//       // cout<<"           Pos trk AODV0 label v0: "<<myvertexv0->GetDaughter(0)<<endl;
//       // cout<<"           Neg trk AODV0 label v0: "<<myvertexv0->GetDaughter(1)<<endl;

//       //   new((fV0s())[nV0++])AliAODv0(* myvertexv0);

//       // AliAODv0* myvertexv0 = new((fV0s())[nV0++]) AliAODv0(vV0, 
//       //  							   rDca3LHDaughters,
//       //  							   rDca3LHToPrimVertex,
//       //  							   rMomPos1,
//       //  							   rMomNeg1,
//       // 							   rDCA);
      
//       new((*fVZERO)[ntracksD++]) AliAODTrack(*HeTrack);
//       new((*fVZERO)[ntracksD++]) AliAODTrack(*PionTrk);

//       // cout<<"           Pos idx from AODV0: "<<myvertex->GetPosID()<<endl;
//       // cout<<"           Neg idx from AODV0: "<<myvertex->GetNegID()<<endl;
//       // cout<<"           Pos trk AODV0 label: "<<myvertex->GetDaughter(0)<<endl;
//       // cout<<"           Neg trk AODV0 label: "<<myvertex->GetDaughter(1)<<endl;
      
//        AliAODv0 *v0test = (AliAODv0 *)fTracks->At(nV0);
//        //cout<<"V0 dal tclones array "<<v0test->GetPosID()<<endl;

//       // cout<<"           Pos idx from AODV0 v0:  "<<myvertexv0->GetPosID()<<endl;
//       // cout<<"           Neg idx from AODV0 v0:  "<<myvertexv0->GetNegID()<<endl;
//       // cout<<"           Pos trk AODV0 label v0: "<<myvertexv0->GetDaughter(0)<<endl;
//       // cout<<"           Neg trk AODV0 label v0: "<<myvertexv0->GetDaughter(1)<<endl;

//       //  delete vertex;
//       //delete myvertex;	
//       //  delete vV0;

//       //cout<<"Fatto vertice secondario AOD"<<endl;
//       //  return;
//     }//pion loop

//     // cout<<"    Out loop Pos from V0: "<<posFromV0<<" He Tracks: "<<HeTrack<<endl;
//     // if(HeTrack && posFromV0!=0)
//     //   {     //     fAODTrackRefs->AddAt(HeTrack,posFromV0);
//     // 	cout<<"if helio: "<<posFromV0<<endl;
//     // 	cout<<"He trk "<<HeTrack<<endl;
//     // 	cout<<"********************************************fAODTrackRefs "<< fAODTrackRefs<<endl;

//     // 	//	HeTrack = static_cast<AliAODTrack*>(fAODTrackRefs->At(posFromV0)); 
//     //   }
//   }

 
//   //  new((*fTracks)[nV0++]) AliAODv0(*myvertex);
//   // delete PionTrk;
//   //delete HeTrack;


//   // lDcaPosToPrimVertex = v0->DcaPosToPrimVertex();
//   // lDcaNegToPrimVertex = v0->DcaNegToPrimVertex();
//   // cout<<"Salvato vertice secondario ESD"<<endl;
//   // delete vertex;
//   // delete myvertex;

//   cout<<"Numero dei v0 prodotti :    "<<  nV0<<endl;
//   cout<<"Numero dei figli prodotti : "<<  ntracksD<<endl;
  
// }




//--------------------------------------------

void AliAODMuonReplicator0::ReconstructAndSaveSecondaryVertex(const AliAODEvent& source, TArrayI &Pions, TArrayI &He)

{

  Int_t ntracks(0);

  AliAODVertex *vtx = source.GetPrimaryVertex();
  AliESDv0 *vertex ;

  //! A set of very loose parameters for cuts 

  Double_t fgChi2max=33.;     //! max chi2
  Double_t fgDNmin=0.05;      //! min imp parameter for the 1st daughter = 500um
  Double_t fgDPmin=0.05;      //! min imp parameter for the 2nd daughter = 500um
  Double_t fgDCAmax=1.;       //! max DCA between the daughter tracks in cm
  Double_t fgCPAmin=0.9;      //! min cosine of V0's pointing angle
  //  Double_t fgRmin=0.2;    //! min radius of the fiducial volume //original
  Double_t fgRmin=0.1;        //! min radius of the fiducial volume = 1 mm 
  Double_t fgRmax=200.;       //! max radius of the fiducial volume = 2 m

  //------------------------------------------
  
  Double_t dcaPionToPV[2] = {1e6,1e6};
  Double_t covarPionToPV[3] = {1e6,1e6,1e6};
    
  Double_t dcaHeToPV[2] = {1e6,1e6};
  Double_t covarHeToPV[3] = {1e6,1e6,1e6};

  Double_t DcaPionToPrimVertex=0;
  Double_t DcaHeToPrimVertex=0;

  //---------------------------------------------------------------------------------------------

  AliAODTrack  *PionTrk = 0x0;
  AliAODTrack  *HeTrack = 0x0;
  
  AliAODv0 *myvertex;

  Int_t nHe    = He.GetSize();
  Int_t nPions = Pions.GetSize();

  cout<<"nHe: "<<nHe<<" npions: "<<nPions<<endl;

  for (Int_t i=0; i<nHe; i++){                            //! He Tracks Loop

    Int_t He3idx=He[i];
    
    HeTrack=source.GetTrack(He3idx);
     
    if (HeTrack){
      HeTrack->PropagateToDCA(source.GetPrimaryVertex(),source.GetMagneticField(),100.,dcaHeToPV,covarHeToPV);
      DcaHeToPrimVertex = TMath::Sqrt(dcaHeToPV[0]*dcaHeToPV[0]+dcaHeToPV[1]*dcaHeToPV[1]);
    }
    
    //-----

    for (Int_t k=0; k < nPions; k++) {                           //! Pion Tracks Loop

      Int_t PionIdx=Pions[k];
      
      PionTrk=source.GetTrack(PionIdx);
      
      if (PionTrk) {
	PionTrk->PropagateToDCA(source.GetPrimaryVertex(),source.GetMagneticField(),100.,dcaPionToPV,covarPionToPV);
	DcaPionToPrimVertex = TMath::Sqrt(dcaPionToPV[0]*dcaPionToPV[0]+dcaPionToPV[1]*dcaPionToPV[1]);	
      }
      
      if ( DcaPionToPrimVertex < fgDNmin)                //OK
	if ( DcaHeToPrimVertex < fgDNmin) continue;      //OK
      
      Double_t xn, xp;
      Double_t dca=0.;
      
      AliExternalTrackParam trackInHe3TPC,trackInNegPion;

      trackInHe3TPC.CopyFromVTrack(HeTrack);  //copy constructor for AliExternalTrackParam object
      trackInNegPion.CopyFromVTrack(PionTrk);

      trackInHe3TPC.CopyFromVTrack(static_cast<AliAODTrack*>(HeTrack));
      trackInNegPion.CopyFromVTrack(static_cast<AliAODTrack*>(PionTrk));	
      
      dca= trackInNegPion.GetDCA(&trackInHe3TPC,source.GetMagneticField(),xn,xp);  //!dca tracks
      
      cout<<"DCA "<<dca<<endl;

      if (dca > fgDCAmax) continue;
      if ((xn+xp) > 2*fgRmax) continue;
      if ((xn+xp) < 2*fgRmin) continue;
      
      //From AliV0Vertex
      
      Bool_t corrected=kFALSE;
      if ((trackInNegPion.GetX() > 3.) && (xn < 3.)) {
	//correct for the beam pipe material
	corrected=kTRUE;
      }
      if ((trackInHe3TPC.GetX() > 3.) && (xp < 3.)) {
	//correct for the beam pipe material
	corrected=kTRUE;
      }
      
      if (corrected) {
	dca=trackInNegPion.GetDCA(&trackInHe3TPC,source.GetMagneticField(),xn,xp);
	if (dca > fgDCAmax) continue;
	if ((xn+xp) > 2*fgRmax) continue;
	if ((xn+xp) < 2*fgRmin) continue;
      }
 
      //===================//
      // "V0" with Tracks  //
      //===================//

      trackInNegPion.PropagateTo(xn,source.GetMagneticField()); 
      trackInHe3TPC.PropagateTo(xp,source.GetMagneticField());
      
      vertex = new AliESDv0(trackInNegPion,PionIdx,trackInHe3TPC,He3idx);
      vertex->SetDcaV0Daughters(dca);
      // vertex->SetDcaPosToPV(DcaHeToPrimVertex);
      // vertex->SetDcaNegToPV(DcaPionToPrimVertex);
      vertex->SetV0CosineOfPointingAngle(vertex->GetV0CosineOfPointingAngle(vtx->GetX(),vtx->GetY(),vtx->GetZ()));

      cout<<"Fatto vertice secondario ESD"<<endl;
      if(vertex->GetV0CosineOfPointingAngle(vtx->GetX(),vtx->GetY(),vtx->GetZ())>fgCPAmin)continue;

      //      cout<<"\n\nVetexer : DCA pos to pv: "<< DcaHeToPrimVertex<<" to be compared to "<<vertex->GetDcaPosToPV()<<" neg "<< vertex->GetDcaNegToPV()<<" cpa: "<<vertex->GetV0CosineOfPointingAngle(vtx->GetX(),vtx->GetY(),vtx->GetZ())<<endl<<endl;

      //return vertex;

      Double_t rDca3LHDaughters = dca;
      Double_t rDca3LHToPrimVertex = vertex->GetD(vtx->GetX(),vtx->GetY(),vtx->GetZ());
  
      Double_t rMomPos[3]={0.,0.,0.};
      Double_t rMomNeg[3]={0.,0.,0.};
      
      Double_t rDCA[2]={0.,0.}; //to be fixed
  
      //ramona 
      
      // rDCA[0] =  vertex->GetDcaPosToPV();
      // rDCA[1] =  vertex->GetDcaNegToPV();
      rDCA[0] =  DcaHeToPrimVertex;
      rDCA[1] =  DcaPionToPrimVertex;
      
      //to be understand if pos is 3He and neg is pi in each case
  
      vertex->GetPPxPyPz(rMomPos[0], rMomPos[1], rMomPos[2]); //is momentum at SV!!
      vertex->GetNPxPyPz(rMomNeg[0], rMomNeg[1], rMomNeg[2]);
  
      const Double_t rMomPos1[3]={2*rMomPos[0], 2*rMomPos[1], 2*rMomPos[2]};
      const Double_t rMomNeg1[3]={rMomNeg[0], rMomNeg[1], rMomNeg[2]};

      // //      AliAODv0 *myvertex=new AliAODv0(vtx,rDca3LHDaughters,rDca3LHToPrimVertex,rMomPos1,rMomNeg1,rDCA);
   
      myvertex=new AliAODv0(vtx,rDca3LHDaughters,rDca3LHToPrimVertex,rMomPos1,rMomNeg1,rDCA);

      cout<<"cpa: "<<vertex->GetV0CosineOfPointingAngle(vtx->GetX(),vtx->GetY(),vtx->GetZ())<<endl;
      cout<<"momPospx: "<<rMomPos[0]<<endl;
      cout<<"momNegpx: "<<rMomNeg[0]<<endl;

      new((*fTracks)[ntracks++]) AliAODv0(*myvertex);
      
      delete vertex;
      delete myvertex;	

      //cout<<"Fatto vertice secondario AOD"<<endl;
    }
  }
  
  // lDcaPosToPrimVertex = v0->DcaPosToPrimVertex();
  // lDcaNegToPrimVertex = v0->DcaNegToPrimVertex();
  // cout<<"Salvato vertice secondario ESD"<<endl;
  // delete vertex;
  // delete myvertex;


}




//______________________________________________________________________________
TClonesArray& AliAODMuonReplicator0::fV0s()
{
  TClonesArray pippo("AliAODv0",30);//= new TClonesArray("AliAODv0",30);
  return pippo;
  //  return *(AODEvent()->GetV0s());
  //  return *(source->GetV0s());
  // return *0x0;
}
